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1. Introduction 

Cepheid and RR Lyrae modelling has a long history going back to the early 
1960's (recently reviewed in Gautschy & Saio 1995). Right from the beginning 
it was quite clear that convection had to be present in the pulsating envelopes. 
Furthermore, modelling showed that convection was necessary to provide a red 
edge to the instability strip, i.e. to stabilize the stars at lower temperatures. 
However convection was deemed to have a minor effect on the shape of the light 
curves and radial velocity curves. And, indeed, purely radiative models gave 
good agreement with the observations of the Galactic Cepheids (e.g. Moskalik 
et al. 1992), although a few problems of varying degree of severity persisted 
(cf. Buchler 1998), such as the inability of radiative codes to model beat pul- 
sations in either Cepheids or RR Lyrae. The light curves of the so-called Beat 
Cepheids or RR Lyrae indicate that these stars pulsate with two basic frequen- 
cies, and with constant power in these frequencies. In addition, radiative codes 
give pulsation amplitudes that are much too large when compared to the obser- 
vations. Furthermore the amplitudes depend on the fineness of the numerical 
mesh. The amplitudes as well as the stability of the limit cycles also depend on 
the values chosen for the pseudo-viscosity. 

In the last few years a wealth of data on variable stars in the Magellanic 
Clouds (MC) has been obtained as a by-product of the EROS and MACHO 
microlensing projects. Because these galaxies have a metal content that is only 
one quarter to one half that of our Galaxy, our observational data base has 
therefore been substantially broadened. Calculations with radiative codes show 
rather clearly that purely radiative models are incapable of agreement with 
observations (e.g. Buchler 1998). 

The fact that resonances among the vibrational modes give rise to observ- 
able effects (e.g. Buchler 1993) can be exploited to put constraints on the pul- 
sational models and on the mass-luminosity relations. The best known of these 
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Figure 1. Time dependence of turbulent energy during the pulsation cy- 
cle. Left: vs. zone index, right: vs. Lagrangean radius. 

resonances occurs in the fundamental Cepheids (Po/P2=2) in the vicinity of a 
period Po=10 days and is at the origin of the well known Hertzsprung progres- 
sion of the bump Cepheids. MC observations (Beaulieu et al. 1995, Beaulieu & 
Sasselov 1997, Welch et al. 1997) show that the resonance may be slightly shifted 
to half a day or a day higher in period. Structure also appears in the Fourier 
decomposition coefficients of the first overtone Cepheid light curves, and is most 
likely due to a resonance P\/P±=2 with the fourth overtone as first pointed out 
by Antonello et al. (1980). Again MC observations indicate that the resonance 
center occurs approximately at the same period. When used to constrain purely 
radiative models (Buchler, Kollath, Beaulieu & Goupil 1996) one obtains stel- 
lar masses that are much too small to be in agreement with stellar evolution 
calculations. 

Improvements to the radiative Lagrangean codes have been made in recent 
years: Adaptive mesh techniques have been used to resolve sharp spatial features 
such as shocks and ionization fronts. Instead of treating radiation in an equi- 
librium diffusion approximation, the equations of radiation hydrodynamics have 
been implemented. However, all these changes have not substantially improved 
the agreement between modelling and observations. It has become patently clear 
that some form of convective transport and of turbulent dissipation is needed if 
we want to make progress. 

Turbulence and convection are inherently 3D phenomena. While a great 
deal of progress has been made in 3D simulations, it remains very difficult to 
model astrophysically realistic conditions which have very large Rayleigh num- 
bers Raw 10 12 , and very small Prandtl numbers Pr~ 10 6 . It is of course even 
more difficult to incorporate them in stellar models (see however the solar models 
of Nordlund & Stein at this meeting) . 

Large amplitude stellar pulsations increase the difficulties by involving time 
dependence. Indeed, the source regions occur in the partial ionization regions 
of hydrogen, helium and Fe-group atoms, and these features are neither La- 
grangean (they move through the fluid) nor Eulerian (they move through space) . 
Fig. [I] shows the behavior of the turbulent energy et over a period in a pulsating 
Cepheid model with a period of 10.9 days (M=6.1M , L=3377L , T e// =5207K, 
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Figure 2. Convective luminosity; left: vs. zone index, right: vs. La- 
grangean radius. 

X=0.70, Z=0.02). Similar behavior occurs in RR Lyrae models, but we concen- 
trate here on Cepheids which are actually more daunting numerically because 
of the sharpness of their H ionization front. On the right side we display et as a 
function of Lagrangean radius. The lightness of the grey reflects the strength of 
et- On the left, Fig. |l| displays et as a function of zone index, i.e. as attached to 
the Lagrangean mass coordinate, i.e. in the fluid frame. The turbulent energy 
is largest in the region associated with the combined H and first He ionization 
fronts. The next most important region of turbulent energy is the He + -He ++ . 
There can also be turbulent energy in the Fe group partial ionization regions, at 
least for Galactic metallicity, but is comparatively weak and does not show up 
on the scale of the figure. Fig. [j] clearly shows how the turbulent energy tracks 
the source regions which move through the fluid during the pulsation. It also 
shows the importance of time dependence in the convective pulsating envelope. 
Both figures show that the turbulent energy increases during the pulsational 
compression phase and that the two turbulent zones briefly merge. 

Fig. |2] shows the temporal behavior of the convective flux in the frame of 
the zone index (Lagrangean) and in the stellar frame, respectively, and can be 
compared to the turbulent energy in Fig. |l|. The convective flux exists only in the 
regions of negative entropy gradient regions (Y > 0, cf. Eq. |||) and is therefore 
confined to narrower regions than the turbulent energy which can diffuse outside 
these regions. 

Attempts at including convection in pulsation codes are actually not new. 
Castor (1971) was the first to present a nonlocal time dependent formulation and 
numerical application to a pulsating stellar envelope model. Later, Stellingwerf 
simplified Castor's formulation and performed a number of model calculations 
(Stellingwerf 1982, Bono & Stellingwerf 1994). Similar turbulent convective 
model equations have been used by Gehmeyr and Winkler (1992). Another 
early computation of linear convective models is that of Gonczi & Osaki (1980). 
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2. The Turbulent Convection Recipe 
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where p is the gas pressure, (3 is the thermal expansion coefficient, A = a\ H p , 
and Hp = pr 2 /(pGM) is the pressure scale height, and other symbols have their 
usual meanings. 

This scheme gives rise to an unphysical behavior at the boundaries of the 
convective regions where Y — » 0. Because 5Sk oc 5Y/y/Y the linearization has 
a pole there that shows up in the growth rates along sequences of models when 
the zoning is very fine or the mesh happens to fall on a point with small Y. 

This difficulty can be avoided with an alternative model equation for the 
source which is more in line with the Gehmeyr- Winkler formulation (1992) 



St 



(a s a A ) 2 V -(3TY , 
P 



(7) 



For comparison, in Fig. ^, we show the effect that the two formulations have on 
the linear growth rates for a sequence of Cepheid models (M=6.75M0, L=4843, 
X=0.70, Z=0.02, variable T e ff). The difference is seen not to be substantial, 
although the alternate St increases the maximum period that unstable overtone 
models can have (cf. Fig. 12 and 13 in YKB). For this model we have taken the 
parameters (a a = 1.0, a c = 2.25, a s = 0.75, a v = 1.8, at = 0.25, a p = 0.667, eo = 
l.e4, a a = 0.4). These values are used for illustrative purposes. Except for this 
one example we have not yet explored the general effects of using expression ([/]) . 
We are also still in the process of calibrating the a's using the best available 
astrophysical constraints. 



3. Work Integrand 

It is interesting to see how turbulent convection affects the stability of the pul- 
sational modes. Turbulent convection actually affects the stability in two ways, 
indirectly, by altering the structure of the equilibrium model and, directly, in 
the linearization of the equations, (cf. e.g. Yecko, Kollath & Buchler 1998, YKB 
hereafter). 
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Figure 3. Fundamental and first overtone growth rates for a Cepheid se- 
quence. Solid line: formulation of Eq. 0, dotted line: Eq. m 



The work done on the pulsation per cycle is given by 

W = jdt(J dmp^), (8) 

where the total pressure p = p g +Pt+Pu is composed of the separate contributions 
of the gas (and radiation) pressure p g , the turbulent pressure p t and the eddy 
viscous pressure p u . If we denote the linear eigenvalues by a = ito + k, for an 
assumed exp(crt) dependence, then the relative growth rate is given by 

k 2tt f 
r\ = 2— = — — Im / dpov am , (9) 



LO UJ 2 I 

1 = J \5r\ 2 dm, (10) 

and the 5 refer to the pressure, specific volume and radial displacement parts of 
the modal eigenvector, respectively, and / is the moment of inertia of the mode. 

The quantity r\ represents the energy growth of a mode over one period, 
equal to the inverse of the quality factor Q that is commonly associated with 
resonant electronic devices. 

Here we illustrate with a fundamental Cepheid model (with M=5.2M , 
L=3293L Q , T e jf=5677K, X=0.716, Z=0.01) how the work integrands are af- 
fected by convection. It is of interest to see how the various regions contribute 
to the work, as well as how the turbulent convective quantities affect the stabil- 
ity. 

In Fig. |I| we display the linear work integrand (thick solid line) together with 
the separate contributions: gas pressure (thin solid line), pt (dotted line) and 
p u (dashed line). The area under the curve is slightly positive since the mode 
is linearly unstable. As expected, the eddy pressure is everywhere damping. 
The turbulent pressure, on the other hand, can be both driving or damping 
depending on its phase with respect to the density variations. The sharp peak 
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Figure 4. Linear Work Integrand (surface to the right) showing separate 
contributions of the gas, turbulent and eddy viscosity pressures. 

is associated with the H and first He ionization, whereas the broad peak is due 
to the second He ionization and the tiny peak to Fe. 

The total work (Eq. that is done over a limit cycle is zero, but again it is of 
interest to see how nonlinear effects change the nonlinear work integrand which 
is displayed in Fig. |5[ Our nonlinear work integrand is arbitrarily normalized 
by twice the nonlinear pulsational kinetic energy. The figure shows the separate 
contributions of p g , pt and p v . Here there is, in addition, a pseudo-viscous 
pressure whose contribution is very small compared to that of the other pressures 
and is not shown here. 

In comparison to the linear work integrand, most noticeable are (a) the 
broadening of the driving region because the (non Lagrangean) ionization fronts 
sweep through the envelope during the pulsation; this is already known from 
purely radiative models (e.g. Figs. 4 and 7 in Buchler 1990) and (b) the greatly 
enhanced damping by the eddy viscosity pressure p u . 

A final comment concerning frequently made approximations. In many early 
pulsation computations convection was assumed to be 'frozen in': Convection 
was included in the computation of the equilibrium model, but all convective 
quantities were held constant in the calculation of the period and of the linear 
growth rates. A similar approximation is often made in stellar evolution compu- 
tations. In YKB we have examined this approximation and found it to be very 
lacking. The perturbation of the turbulent quantities, and concomitantly of the 
convective flux, has a very strong damping effect on the pulsation. 

In Fig. | we summarize these results for a resequence of models (with 
M=5Mq, L=2060L q ) for the fundamental and first overtone modes. The solid 
lines represent the exact growth rates (i.e. correct linearization of all quantities). 
The line with crosses represents the 'frozen convection' approximation which is 
seen to be inadequate. The fundamental instability strip (the domain where 
the modal growth rate is positive) is enormously broadened and shifted. For 
the overtone the effect appears even more drastic. The mode, which is stable 
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Figure 5. Nonlinear Work Integrand (surface to the right) showing sepa- 
rate contributions of the gas, turbulent and eddy viscosity pressures. 

throughout the whole temperature region, now becomes unstable over a very 
broad region. 

The dotted line corresponds to the approximation of 'adiabatic' convection, 
i.e. T5st = Set — pt$v = 0. Physically it corresponds to assuming that all con- 
vective time scales are very long. This approximation is seen to underestimate 
somewhat the damping effects of convection. 

Another convenient approximation is the other extreme, which is to assume 
that all convective time scales are very short compared to the other time scales, 
i.e. from Eq. || we obtain 

£ (r 2 F t ) - e \a d (e t - St) = 0. , (11) 

This is seen to be the best of the approximations. It is also the simplest to apply 
in evolutionary calculations in which a time independent, local mixing length 
recipe is used (which would correspond here to setting in addition Ft = or 
at = 0, i.e. no diffusion of turbulent energy). 

4. Nusselt vs. Rayleigh Numbers 

The Nusselt number is defined as Nu=F c / F con d, where in our case the conductive 
flux is the radiative flux, and the Rayleigh number is ~R&=gPd 3 TY/(i/x)- Here 
g is the local gravity, d is the local scale height, v is the kinematic viscosity 
and x is the radiative conductivity. There is general agreement that Nu should 
depend on Ra, viz. Nu = Ra a , but there is no theoretical agreement on what 
the value of a should be (e.g. Spiegel 1971). Some experimental results indicate 
that a = 0.28 (Castaing et al. 1989), but it is not clear that they should apply to 
the stellar case where the boundaries can adjust to accommodate a fixed stellar 
luminosity, and where the physical quantities have strong spatial variations, 
especially through the partial ionization zones. 
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Figure 7. Left: Local Nusselt versus Rayleigh number. Right: Peclet 
Number in two Ccpheid envelopes. 
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In Fig. [?] we reproduce that behavior of the local Nu versus Ra numbers 
throughout the convective regions of the two typical Cepheid models M=5M Q , 
L=2O9OL , T e// =4900 (solid line) and 5300K (dotted) of YKB. Only the com- 
bined H-He convective regions, where Nu> 1, are shown. For reference we have 
shown two thin lines with slopes 1/2 and 1/3, respectively. Throughout the 
convective region the exponent a varies between 0.45 and 0.53, and thus agrees 
best with the higher theoretical value of 1/2 (Spiegel 1971). The right hand side 
shows the Peclet number defined as the ratio of the thermal diffusion time scale 
to the convective time scale. 

5. Sequence of Cepheid models 

YKB performed some sensitivity tests of the properties of Cepheid models ob- 
tained with the ID turbulent convective model. Here we just add Fig. ^ which 
displays the strength of the convective luminosity as a function of zone number 
(bottom scale), for a sequence of Cepheid models starting from the blue edge 
in front, to the red edge in the back. The importance of the convective flux 
increases from the blue edge, where it is relatively unimportant, to the red edge. 
The H and first He ionization regions are always merged into a single zone. Near 
the blue edge the second ionization region for He forms a separate convective 
region, but when we arrive at the red edge, convection encompasses both H and 
He regions, and almost joins with Fe region (left). 

6. Double-Mode (DM) Pulsations 

The numerical modelling of double mode (DM) pulsations has been a long stand- 
ing quest in which purely radiative models have failed. In a recent paper (Kollath 
et al. 1998, hereafter KBBY) it was shown that with the inclusion of turbulent 
convection DM pulsations appear almost naturally in Cepheid models. Almost 
concomitantly, but independently, Feuchtinger (1998) found DM behavior in RR 
Lyrae pulsations which we have since also confirmed. KBBY described the be- 
havior of the DM Cepheids in terms of truncated amplitude equations (Eqs. 1 
of KBBY), and they appeared to give excellent agreement with the model that 
was studied. 

Fig. 1 of KBBY showed the transient evolutions for a given Cepheid model 
and for different initializations of the hydrocode. The evolution toward a DM 
pulsational state is clearly exhibited. The results of the pulsational states of a 
number of Cepheid models were summarized in a bifurcation diagram (Fig. 4 
of KBBY) . The DM states were obtained with the regular hydrodynamics code 
after lengthy time integrations with suitable initial conditions. The single mode 
pulsational states, whether stable or not, were obtained with Stellingwerf's relax- 
ation method (cf. Kovacs & Buchler 1987), sometimes with a lot of perseverance. 

When such models were more carefully scrutinized it became apparent that 
a different transient evolution was possible (Kollath et al. 1999), namely toward 
the F limit cycle on the bottom right. This situation is shown in Fig. ||. It is 
clear that in addition to the stable DM there must coexist a stable F limit cycle 
and a second unstable DM. 
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Figure 8. Profiles of the ratio of convective to total luminosity along a 
sequence; blue edge in front (bottom) and red edge towards the back (top). 
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Figure 9. Transient evolution of Cepheid model. Left (right): the F am- 
plitude - 01 amplitude plane; right: energy plane The trajectories correspond 
to various initializations. 
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Figure 10. Bifurcation diagram for the Cepheid model sequence; 
Solid/open circles: F and 01 limit cycles, Squares: F and 01 component 
amplitudes of double mode cycles. 



This new development forces us to also reconsider the bifurcation diagram 
(Fig. 4 of KBBY). We had suggested that the sharp vertical rise (drop) of the 
fundamental (overtone) amplitudes was due to the presence of a pole in the 
discriminant T>. While such a pole is present it turns out to be too far away (in 
T e ff) to cause the observed vertical slope. Upon closer inspection it has been 
found (Kollath et al. 1999) that the bifurcation diagram is a bit more complicated 
than first thought. Fig. 1C is an adaptation of the results of Kollath et al. (1999). 

Indeed, the region of fundamental mode pulsations extends to the left into 
the region where DM pulsations can occur. There is thus a narrow region of 
hysteresis where both F and DM pulsation can occur. We note immediately that 
this bifurcation structure, in particular the hysteresis, cannot be accommodated 
with the amplitude equations of KBBY that were truncated at the A 3 terms. 

Kollath et al. (1999) show that one can readily get agreement by adding 
the most important next order terms in the truncation, which normal form 
theory shows to be —tqAq and —r\A\. (We disregard the additional quintic 
cross-coupling terms). 



dA 



o 



dt 
dAi 

dt 



(«o - QooAl - qoiAj - tqAq) A 
(ki - q w A 2 - q n Aj - r x A\) A x 



(12) 
(13) 



Rather than using amplitudes A, it is equivalent and perhaps more convenient 
here to introduce the 'energies', e = A 2 , instead of the amplitudes A. The 
amplitude equations, with the new terms added, take on the form 



— = 2 (ko - gooeo - 901 ei - r e ) e 

— =2(ki- qi e - q n e l - nej e x 



(14) 
(15) 
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Figure 11. The fixed point solutions as a function of the control pa- 
rameter A. Solid dots for stable, open circles for unstable fixed points. 



The loci of the fixed points are obtained by setting the RHSs of these equations 
equal to zero. Without the r terms these nullclines (other than the two coordi- 
nate axes) are simply straight lines that intersect at most once, and when they 
intersect they give the DM as has been known for a long time (Buchler & Kovacs 
1986). Clearly no hysteresis is possible in this case. On the other hand, even 
with small r values the lines bend and multiple intersections become possible. 



This situation is depicted in Fig. 11. For the sake of illustration, we have 



chosen the numerical values ((?oo=2.179e-3, (?oi = 4.5e-3, q , io=5.9e-3, gn=16.e-3, 
ro=-3.e-4, ri=1.4e-3), for simplicity keeping these values constant even though 
in a real sequence of models they would vary. The variation of the growth rates 
along a sequence is more important and we assume that kq = Rq + 0.8e-3A, K\ = 
R% + 1.6e-3 A, where Ro=2.e-3, Ki=7.5e-03. The parameter A varies between 
and 1 along this sequence. 



The corresponding bifurcation diagram is presented in Fig. 12. It is seen to 



display the same general features as the actual Cepheid diagram. In particular, 
it has an single-mode 01 regime up to A ~ 0.08, a DM regime from A ~ 0.08- 
0.90, and a coexistence between DM and F modes from A ~0. 67-0. 90. To the 
right A ^ 0.90, only the F mode LC is stable. Note that the annihilation of 
the stable and unstable DMs that occurs at A ~ 0.90 gives rise to the vertical 
observed tangent. 

Note that the complexity of the bifurcation diagram is partly due to the 
values we have chosen for the control parameters, T e jf and a u . Our values cor- 
respond to realistic Cepheids and are not idealized for the purpose of clarifying 
the evolution into single and double modes. If we had chosen instead to unravel 
the complete nature of the bifurcation, we would have been forced to choose 
both T e ff and a v to correspond to the polycritical point - where the F, O, DM 
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Figure 12. Bifurcation diagram corresponding to the illustrative ex- 
ample. Amplitudes of the SM and DM limit cycles. The stable (unsta- 
ble) limit cycles are denoted by thick (thin) dots. 



and the trivial solution coexist (this point was previously discussed in KBBY 
and plotted there in Fig. 3). Near to this point the dynamics is given by the 
cubic equations (we can take this as the definition of near). Furthermore, the 
bifurcation structure is straightforward once we know how we move through the 
parameter space given by T e ff and a u . (This is not true if the bifurcation is 
subcritical, but as yet we have not encounterd this case). For a more general 
unfolding of the bifurcation, however, this ideal picture is easily extended be- 
yond its reach. So we ought to expect that some effects of this breakdown, in 
the form of an increasing nonlinearity, will begin to appear. What we seem to 
be witnessing here is, in fact, the need for quintic terms as the polycritical point 
becomes more distant. 



7. Conclusions 

It is perhaps remarkable that such a simple ID recipe for turbulent convection 
can give such drastic improvements over purely radiative codes. It may indicate 
that, at least for Cepheid and RR Lyrae variables, this recipe incorporates all 
the physics of turbulence and convection that is essential to model these pul- 
sations. It is possible that a nonlocal, time dependent dissipation which this 
model equation provides is all that is needed. 

We are only at the beginning of the process of calibrating the seven a 
parameters that appear in the turbulent convective description. There are nu- 
merous constraints that need to be satisfied, and we hope that despite the large 
number of these a's these constraints can be satisfied. In particular it will be a 
challenge to obtain the observational properties of both the Galactic and of the 
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Magellanic Cloud variable stars. Only then will we know whether our simple ID 
model is adequate. 
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